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ABSTRACT 

A multi-dimensional, coupled thermal response modeling system for analysis of hypersonic entry vehicles is presented. 
The system consists of a high fidelity Navier-Stokes equation solver (GIANTS), a two-dimensional implicit thermal 
response, pyrolysis and ablation program (TITAN), and a commercial finite-element thermal and mechanical analysis 
code (MARC). The simulations performed by this integrated system include hypersonic flowfield, fluid and solid 
interaction, ablation, shape change, pyrolysis gas^generation and flow, and thermal response of heatshield and structure. 
The thermal response of the heatshield is simulated using TITAN, and that of the underlying structural is simulated 
using MARC. The ablating heatshield is treated as an outer boundary condition of the structure, and continuity 
conditions of temperature and heat flux are imposed at the interface between TITAN and MARC. Aerothermal 
environments with fluid and solid interaction are predicted by coupling TITAN and GIANTS through surface energy 
balance equations. With this integrated system, the aerothermal environments for an entry' vehicle and the thermal 
response of the entire vehicle can be obtained simultaneously. Representative computations for a flat-faced arc-jet test 
model and a proposed Mars sample return capsule are presented and discussed. 

INTRODUCTION 

Spacecraft heatshields typically use thermal protection system (TPS) materials that pyrolyze and ablate at high 
temperature for mass-efficient rejection of the aerothermal heat load. Pyrolysis is an internal decomposition of the solid 
that releases gaseous species, whereas ablation is a combination of processes that consume heatshield surface material 
(including chemical reactions, melting, and vaporization). For design and sizing of ablating spacecraft TPS materials, 
it is important to have a reliable numerical procedure that can predict surface recession rate, in-depth pyrolysis, and 
internal temperature history. The thermal properties and the boundary conditions of the heatshield also must be 
precisely defined to perform an accurate prediction. 

A simplified schematic diagram of rhe geometry is shown in Figure 1. The outer face of the heatshield is exposed to 
an aerothermal heating environment that can be predicted using a high fidelity computational fluid dynamics (CFD) 
code with surface-thermochemistry' boundary conditions. The inner face of the heatshield is attached to a structure that 
contains various other components The structure often should be included in the thermal response simulation to 
correctly predict the bond-line temperature history'. The maximum allowable bond-line temperature is typically the key 
parameter that drives the TPS thic<ness. The structure and its interior can have a complicated geometry, but the 
thermal response is relatively simple compared with that of the TPS. 

A number of programs have been written to simulate ablating heatshields. The CMA code, 1 developed by Aerotherm 
Corporation in the 1960’s, was one of the first one-dimensional codes. CMA solved internal energy balance and 
decomposition equations coupled with general surface energy balance boundary (SEB) conditions to simulate the 
response of pyrolyzing and ablating heatshields in hypersonic flows. In 1997, the FIAT code was developed at NASA 
Ames Research Center to support the development of lightweight ceramic ablators. FIAT is numerically more stable 
and solves a wider range of problems than CMA, and it has been used for TPS sizing calculations in various NASA 
space missions. 2 However, the governing equations solved in both CMA and FIAT are one-dimensional. Thus, 
neither of these two codes is applicable for conditions in which a one-dimensional assumption is not true, such as the 
nose-tip of a slender hypersonic reentry vehicle or around the shoulder of a blunt vehicle. 

For axisymmetric or planar geometries, Aerotherm also developed ASCC, 3 an engineering code for fully-coupled 
fluid/solid simulation of ablating heatshields. In ASCC the thermal diffusion equation is solved using a finite- 
difference scheme with overlaid grids. A two-dimensional ablation code developed at Sandia National Laboratories 
uses the finite-control-volume method with unstructured grids. 4 In this code, mesh motion (owing to ablation) is 
modeled by assuming the mesh behaves as a linear elastic solid. In both of these codes the effect of pyrolysis gas is 
not considered. Recently the THAN program was developed to perform high fidelity two-dimensional thermal 
response and shape change simulations for pyrolyzing ablators. 5 The governing equations, which include energy 
conservation and a three-component decomposition model with a moving grid, are discretized using a finite-volume 
approximation and general body-fitted coordinates. Time-accurate solutions are achieved by an implicit time marching 
technique using Gauss-Siedel line relaxation with alternating sweeps. 6 In many cases, a coupled fluid and solid 
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simulation is required for accurate prediction of shape change, because shape change can not be correctly determined 
without appropriate aerothermal heating, but aerothermal heating is very' sensitive to the geometry of the solid. To this 
end, the Navier-Stokes solver, GIANTS, 7 was successfully coupled with TITAN using a loosely coupled method for 
simulating the fluid sol id interaction. TITAN uses an external zone with a moving grid and internal grid that is fixed. 
This two-zone approach can hand e some simple geometries but is not sufficiently flexible to model a complicated 
interior structure. 

To address this deficiency, the purpose of the present work is to integrate GIANTS and TITAN with the MARC finite- 
element code. 8 The three codes are integrated through the use of MARCs user-supplied boundary condition 
subroutines. TITAN is used for the heatshield analysis because the physics of pyrolysis and ablation is too 
complicated for available commercial finite-element codes to accurately simulate. The structural materials beneath the 
heatshield experience relatively benign heat loads, and are usually multi-dimensional with complicated geometry. A 
commercial finite-element code, such as MARC, can perform thermal and structural analyses of the interior of a space 
vehicle. This integrated system can predict, simultaneously, two-dimensional aerothermal heating, ablation, shape 
change, and thermal response of heatshield and structure. The integrated simulation program should significantly 
reduce the computational and human effort requfred to design and analyze a spacecraft heatshield and attached structure. 
To demonstrate the capability of this modeling system, representative computations for an arc -jet model and a proposed 
Mars sample return capsule are presented and discussed. 

INTEGRATED MODELING SYSTEM 

Figure 2 presents the flow chart for integrated thermal response modeling system. The MARC, TITAN, and GIANTS 
codes are the three major comporents in the system. A brief description of each code is listed in the following 
subsections. The accuracy of each code has been studied previously. 5 7 9 This modeling system is integrated around the 
MARC code. The front end of MARC, known as Mentat, 10 is a graphical user interface program. The material 
properties are provided by TPSX, 11 a material properties database developed at NASA Ames Research Center. 

MARC and TITAN are interfaced hrough MARCs heat-flux boundary' condition routine. The continuity' conditions of 
heat flux and temperature are imposed at the interface: 

Qmarc = _ kVT | titan (0 

= t marc ( 2 ) 

At each time increment (from / to / + 1) the front-face heat flux boundary condition of MARC is updated using the 
temperature gradient calculated by TITAN using Eq. (1), and the back-face temperature boundary condition of TITAN is 
equal to the temperature computed by MARC using Eq. (2) at each interface point. 

TITAN and the flow environment code GIANTS are interfaced using MARC's surface energy balance boundary 
condition. The conditions at the ablating surface are determined by convective and radiative heating and by surface 
thermochemical interactions with the boundary' layer gases. The SEB equation employed is of the convective transfer 
coefficient type and takes the form: 2 

Pe u e C h[ H r~( > + B ')K ] f m cK + "'A + <?/?*■ ~ Foe w T w ~ <7 CV = 0 (3) 

The sum of the first three terms in Eq. (3) is the convective heat flux plus the heat of ablation. The fourth and fifth 

terms are the radiative heat fluxes absorbed and reradiated by the wall, respectively, and the last term represents the rate 
of conduction into material. Here 8' is a nondimensionalized mass blowing rate. Tables of B' for charring materials 
can be generated using ACE 13 or MAT. 14 

A blowing correction accounts for the reduction in transfer coefficients due to the transpiration or blowing effect of 
gases from pyrolysis and surface ablation being injected into the boundary layer. The blowing correction equation for 
convective heat transfer is: 

C h _ln(l + 2Afl') 

Cm " 2XB' ( ) 

where A is blowing reduction parameter and C h jC h \ is the ratio of blown (ablating) to unblown (nonablating) heat 
transfer coefficients. With A = 0.5 Eq.(4) reduces to the classical laminar-flow blowing correction. 15 A variable A is 
used for cases with transitional or turbulent flow. 

The computation starts with the MARC code. At each time increment, MARC calls TITAN to compute a new heat 
flux boundary condition. At the same time TITAN obtains its back-face temperature boundary condition from MARC 
and, if necessary, calls the flow environment code to update the front-face SEB boundary' condition. When the 
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maximum local surface recession exceeds a predefined criterion since the last surface convective heating was calculated, 
a new' flow-field grid is generated based on the current body geometry, and the flow' simulation routine is called to 
compute the aerothermal heating environment. Each call to the flow' environment routine is a steady-state calculation. 
The non-ablating surface heating is calculated by the flow code, and the blow ing reduction parameter is used in TITAN 
to take into account the laminar flow blockage due to surface mass blowing. As expected, the GIANTS calculation is 
much more computationally intensive than the thermal response computation. For a coupled GIANTS/TITAN/MARC 
simulation, most of the CPU time is consumed by flow environment calculations. Thus, the CPU time required for a 
simulation is primarily determined by the efficiency of the GIANTS code. 

MARC 

Structural material thermal response and thermal stress analysis are performed using the commercial finite-element code, 
MARC. The heat transfer and mechanical analyses can be coupled or performed separately. The internal energy 
conservation law' is 

pc p ~~ - v ‘ (kVT) - Q . * (5) 

' dt 

The boundary conditions include prescribed temperature, heat flux, and convective heat transfer coefficient to the 
environment. For transient analysis, the initial temperature is specified to start the calculation. To perform thermal 
stress analysis, MARC uses an instantaneous thermal expansion coefficient defined as 

d £ ki = a u dT . ( 6 ) 

In most cases, however, thermal expansion data are given w ith respect to a reference temperature 7 q, as 

E ,h -a(T-T 0 ). (7) 

Hence, the conversion of the expansion data to the instantaneous thermal expansion coefficient becomes 
— da _ 

«=«+— (T-T 0 ). (8) 


TITAN 


Ablating heatshield thermal response and shape change computations are performed using TITAN. The governing 
equations include energy conserva ion and three-component decomposition model. The SEB condition, Eq. (3), is 
solved with a moving grid to calculate the shape change due to surface recession. The internal energy balance is a 
transient thermal conduction equation with additional pyrolysis terms: 5 

cfT 

pc p — - v • (k VT) -(hg -h )V* m g + m g • V h g + spc p VT (9) 


The individual terms in Eq. (9) may be interpreted as follows: rate of storage of sensible energy, net rate of thermal 
conductive heat flux, pyrolysis energy-consumption rate, net rate of energy convected by pyrolysis, and convection rate 
of sensible energy due to coordinate system movement. 

A three-component decomposition model is used. The resin filler is presumed to consist of two components which 
decompose separately, while the reinforcing material is the third component which can decompose. The instantaneous 
density of the composite is given by: 

P -npA +P*) + U -T)p c (10) 


where A and B represent components of the resin, and C represents the reinforcing material. T is the volume fraction of 
resin and is an input quantity. Each of the three components can decompose following the relation: 
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where p vi and p ri and are the original and the residual (or terminal) density, respectively, of component /. The motion 
of pyrolysis gas is assumed to be one-dimensional along the surface-normal (q) direction, and thus the mass flow rate 
of pyrolysis gas at the surface is calculated using the following approximation: 
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The governing equations are discretized using a finite-volume approximation with a general body fitted coordinate 
system. A time accurate solution is achieved by the implicit time marching technique using Gauss-Siede! line 
relaxation with alternating sweeps. The computational grid is compressed during the course of computation to account 
for surface recession. 


GIANTS 


The flow simulation over a large angle blunt body is performed using a Navier-Stokes solver. The GIANTS code 
solves the time-dependent conservation equations of mass, momentum, and energy for chemical and thermal non- 
equilibrium flow-field. The species mass conservation equation is given by 6 


dPj d d 

the momentum conservation is written as 

d d dx it. * 

—(pu l ) + —(pu l u k ) — , 

at ax f. ax^ 


and the energy conservation as 
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(15) 


The governing equations are discretized using a finite-volume method. The numerical method used to solve these 
discretized equations is exactly the same as that used in the TITAN code. This fully implicit, Gauss-Seidel line 
relaxation technique has been shown to yield steady-state results efficiently. A bifurcation diffusion model 16 was also 
implemented in this version of GIANTS to correctly compute multi-species mass diffusion. 


ARC JET TEST MODEL 

The first test case is a simulation of the flat-faced cylindrical arc-jet test model shown in Figure 3. The TPS material is 
phenolic-impregnated carbon ablator 17 (PICA) with a body radius of 3.81 cm and comer radius of 0.38 1 cm. The PICA 
is bonded to a piece of aluminum honeycomb structure, and this assembly is inserted into a copper model holder. The 
copper is protected from high temperature by a ring of TPS tile material (AETB-12) and disk of TPS blanket material 
(FRSI). The length of arc -jet heat pulse is 35 s. The TPS is sized to keep the structure below a specified temperature 
for several minutes. 

The PICA material is divided into 2 computational zones. Zone I is the region right next to structural materials, where 
temperature should remain sufficiently low so that resin decomposition can be ignored. Zone II is the region directly 
exposed to hot arc-jet stream, where surface recession and pyrolysis gas effects must be considered in the simulation. 
PICA zone II is finite-volume meshed for TITAN; the rest of the model is finite-element meshed for MARC. Because 
a large subsonic region forms in front of the flat-faced cylinder, the flowfield is simulated using the Navier-Stokes 
solver GIANTS. The flow simulation is limited to the region in front of the PICA for simplicity. The exterior 
surfaces of the AETB-12 and copper are assumed to be adiabatic. 

Predicted heat flux distributions (from GIANTS) along the PICA surface at various times are presented in Figure 4, and 
the recession history (from TITAN) is shown in Figure 5. The predicted Mach contours at the beginning and the end of 
the heat pulse are depicted in Figure 6. Initially the heat flux is very high on the relatively sharp comer and 
significantly lower at the centerline. Therefore the ablation rate is highest near the comer, and the comer gradually 
rounds off. As the comer rounds, t ie bow' shock moves closer to the surface, and consequently the heat flux around the 
stagnation region increases. 

Fluid and solid temperature contours at 35 s are shown in Figure 7 (for clarity, different temperature scales are used for 
the fluid and solid). The maximum shock and solid temperature reaches about 6800 K and 2900 K, respectively. At 
200 s there is no flowfield because the arc jet is turned off, so only the solid temperature is plotted. The PICA surface 
continuously cools down through surface radiation, and as heat is conducted into the model the maximum temperature 
is at an interior location. Thus, at 200 s the hottest spot (about 650 K) is located at a depth of about 0.6 cm from the 
front surface. 


MARS SAMPLE RETURN VEHICLE 

The second test case is simulation of a Mars sample return capsule. The purpose of the Mars Sample Return Mission 
is to return samples of material collected at Mars to Earth so that they can be studied here. The final phase of the 
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mission uses an Earth entry, descent, and landing capsule (Figure 9) that has a 0.9-m diameter, spherically blunted, 60° 
half-angle cone forebody. This Earth Entry Vehicle (EEV) is designed to transport the samples safely through the 
Earth’s atmosphere to a recoverable location on the surface, while maintaining the Martian samples at a temperature 
below 323 K. To minimize system complexity and the number of potential failure modes, the EEV design relies on a 
non-parachuted capsule that decelerates solely by aerodynamic forces during entry'. A passive energy-absorbing material 
is used to cushion the Mars samples from excessively large deceleration during landing. The assumed trajectory is 
shown in Figure 8. The entry velocity is 11.6 km/s at angle of -25°. Peak heating occurs at 17 s, and peak 

deceleration is about 130 g at 20 s Ground impact (’’landing'*) occurs at about 41 m/s at 360.5 s. 

A complete Pro/E model of the EEV contained 57 different types of parts including numerous small, non-ax ismmetric 
features such as antennas, bolts, cable connections, latches, springs, and vents. 19 Suppression of these features resulted 
in the model, shown in Figure 9, that contains 19 axisymmetric parts and volumes. The samples are located in region 
19 inside an aluminum sample canister, parts 13-14. This can is held within a titanium containment vessel, parts 15- 

lb. This assembly resides in a spherical structure called the impact shell, parts 4-5. This shell is made of carbon- 

carbon composite, as is the rest the vehicle substructure, parts 1-3. Volume 17 contains an energy absorbing composite 
structure that is assembled in a soccer-ball construction. The remainder of the vehicle interior is filled with various 
insulating materials, air gaps, and silicone rubber. The forebody heatshield is 1.2-cm thick chopped-molded and tape- 
wrapped carbon phenolic, and the aft heatshield is 1.0-cm thick SLA 561-V. Materials properties were obtained from 
NASA Langley and from the TPSX database. 

The TITAN code is used for analysis of the heatshield materials only, and the rest of the vehicle is solved in MARC. 
The MARC material map is presented in Figure 10, and the three grids are illustrated in Figure 11. A very fine grid is 
needed in both the fluid and the solid to resolve the steep temperature gradients and complicated physics near the 
surface of the vehicle. A relatively coarser grid is adequate in the interior of the vehicle where temperature gradients are 
lower and time scales are longer. 

Predicted forebody pressure distributions at selected times are shown in Figure 12. The peak pressure, which occurs at 
about 20 s, is about 1.15 atm. The pressure decreases gradually on the nose, with a slight minimum near the sphere- 
cone tangency point. The pressure increases slightly along the conical frustrum, then decreases rapidly as the flow 
accelerates around the shoulder of the vehicle. Predicted forebody convective heat flux distributions are presented in 
Figure 13. The heat flux decreases steadily, except for a local maximum near the shoulder during part of the trajectory. 
The peak heat convective flux is about 1.1 kW/cirf at 17 s. The stagnation point heat flux is shown in Figure 14. The 
vehicle decelerates rapidly owing to its bluntness, thus the heat pulse is fairly narrow, with 95% of the heating 
occurring in 14 s. The radiative heat flux over the forebody is estimated using the NOVAR code. 20 The peak radiative 
heat flux is about 0.22 kW/cirf at 16.4 s. The total heat load at the nose is about 9.2 kJ/cm 2 of which 8.5 kJ/crrf is 
convective and 0.7 kJ/cm 2 is radiative. 

The forebody heat flux and pressure histories are imposed as boundary conditions of the integrated TITAN/MARC 
code. Owing to the large uncertainties in predicting afterbody heating, afterbody flowfields were not simulated using 
the Navier-Stokes solver. Instead, the heating is assumed to be uniform over the entire aft surface and equal to 5% of 
the time-dependent forebody stagnation point heating. This procedure is imprecise but possibly conservative in 
estimating the overall heat load experienced by the aft heatshield. 

The blowing reduction parameter cf 0.5 is used in the boundary conditions, Eqs. (3-4) for TITAN. Predicted peak 
surface temperature and total recession are shown in Figure 15. The peak temperature is 2500 to 3200 K over most of 
the forebody. However, the total recession is relatively low (under I mm) owing to the high density of the carbon 
phenolic heatshield. Because the surface recession is much smaller than the nose and shoulder radii, the effect of shape 
change on the aerothermal heating is insignificant for this problem. The surface and bond line temperatures in the TPS 
at the stagnation point are presented in Figure 16. The surface temperature increases with total heat flux, but after peak 
heating the temperature decreases more gradually by reradiation as heat is conducted back toward the surface from the 
interior of the TPS. The bond line lemperature begins to rise at 50 s, and has reached 500 K at 300 s. The stagnation 
point blowing rates are shown in Figure 17. The carbon phenolic initially pyrolyzes in response to the aerothermal 
heating, then when the heat flux is sufficiently high filiation also occurs., 

/y ’ ■" F<i c e r-cc* s* < c>\ *=> 

In-depth temperature contours computed by TITAN and MARC at various times are presented in Figure 18. At 40 s, 
all the heat is contained within the heatshield. The peak surface temperature has cooled to 950 K, but the internal 
temperatures are unchanged from the assumed initial value of 250 K. At 100 s, most of the thermal energy is still 
confined within the heatshield, but some heat has penetrated into the forebody structure and foam insulations. At later 
times the TPS continues to cool by reradiation and conduction to the interior. The air gap at the nose is effective in 
preventing heat conduction, and the main conduction paths are along the skirt structure and through the body foam. A 
small amount of heat conduction through the aft TPS is also visible at the later times. At 500 s (if the vehicle is 
intact) the sample canister is safely below the design limit of 323 K. 
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The vehicle impacts the ground at about 360 s. Depending on the angle of the collision and the ground characteristics 
at the impact site, the containment vessel is displaced some distance through the energy absorbing material toward the 
nose structure, and other vehicle damage can occur as well. It is not obvious what assumptions are conservative for 
continuing the thermal analysis ir the post-landing (pre-recovery) phase. Various scenarios can be simulated using 
MARC, including an impact simulation using assumed mechanical properties for the various materials. Future work 
will consider contact resistances between materials, and mechanical analyses such as calculation of thermal and 
deceleration stresses in the carbon-carbon structure. 

CONCLUSIONS 

A multi-dimensional fully couple*.! fluid/heatshield/structure thermal response modeling system has been developed. 
The system consists of a two-dimensional implicit thermal response and ablation program (TITAN), a commercial 
finite-element thermal and mechanical analysis code (MARC), and a high fidelity Navier-Stokes equation solver 
(GIANTS). This system can simultaneously predict aerothermal environments, heatshield ablation and shape change, 
and structural material thermal response for hypersonic entry vehicles as well as test models in hypersonic ground 
facilities. Representative computations for 2fh arc-jet model and a proposed Mars sample return capsule were 
demonstrated. This integrated syslem can substantially reduce the computational effort required for multi-dimensional 
thermal analysis of ablative heatshield and structural materials, and thus should simplify the design cycle for future 
vehicles. 
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NOMENCLATURE 


A - area, nT 

S’ = ml p e it e C h , nondimensional mass flux 
B a = pre-exponential constan: in Eq.(l 1), s' 1 
C h = Stanton number for heat transfer 
c p - heat capacity, J/kg-K 

E = internal energy in Eq.(1.5), J/kg 

E a = activation energy in Eq.f l 1), J/kmol 
F = view factor 

g = Earth standard gravity, 9.8 1 m/s 2 
H r = recovery enthalpy, J/kg 

h_ = enthalpy, J/kg 

h = partial heat of charring in Eq.(9), J/kg 
k = thermal conductivity, W/m-K 

m = mass flux, kg/m 2 -s 

p = pressure, N/m 2 

Q = heat source term, J/m 3 

q - heat flux, W/m 2 

q c = conductive heat flux, W/m 2 

q R - radiative heat flux, W/m 2 

R = universal gas constant, Jkmol-K 
R h - maximum body radius, n 

R c = comer radius, m 

R n - nose radius, m 

j = stream length, m 

s = surface recession rate, m s 

T = temperature, K 

t = time, s 

u = velocity, m/s 


Vj z species diffusion velocity, m/s 

vv mass source term, kg/m' 

x general coordinate variable, m 

a -■ surface absorptance 

a thermal expansion coefficient, 1/K 

e T - surface emissivity 

e th -■ thermal strain 

T volume fraction of resin 

A = : blowing reduction parameter 

p = : total density, kg/m 3 

a - Stefan-Boltzmann constant, W/m : -K J 

r ~ CFD stress, N/m 2 

= : decomposition reaction order 

rj = : surface -normal coordinate, m 

subscripts 
c = char 

e ~ boundary- layer edge 

g = ; pyrolysis gas 

i = density component (A, B, and C) 

j = : gas species 

k, 1 = coordinate directions 

m, n = index of computation cell 

v =• virgin 

w = wall 

0 = reference value 

superscript 

1 = index of time level 
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Figure 6. Mach contours for flow past arc-jet model at beginning (0 s) and end (35 s) of heating. 
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Figure 8. EEV trajectory. 
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Figure 9. Axisymmetric EEV geometry. 
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Figure 11. GIANTS grid (blue), TITAN grid (green), and 
MARC grid (red). 
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Figure 12. Surface pressure distribution at selected times 
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Figure 13. Surface heat flux distribution at selected times. 
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Figure 15. Peak temperature and total recession profiles. 
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Figure 16. TPS temperature histories at stagnation point. 




Figure 14. Stagnation point heat flux. Figure 17. Surface mass flux at stagnation point. 
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Figure 18. EEV temperature contours at selected times, K. 
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